Causal association of menstrual reproductive factors on the risk of osteoarthritis: A univariate and multivariate Mendelian randomization study

Objective Several observational studies have revealed a potential relationship between menstrual reproductive factors (MRF) and osteoarthritis (OA). However, the precise causal relationship remains elusive. This study performed Mendelian randomization (MR) to provide deeper insights into this relationship. Methods Utilizing summary statistics of genome-wide association studies (GWAS), we conducted univariate MR to estimate 2 menstrual factors (Age at menarche, AAM; Age at menopause, AMP) and 5 reproductive factors (Age at first live birth, AFB; Age at last live birth, ALB; Number of live births, NLB; Age first had sexual intercourse, AFSI; Age started oral contraceptive pill, ASOC) on OA (overall OA, OOA; knee OA, KOA and hip OA, HOA). The sample size of MRF ranged from 123846 to 406457, and the OA sample size range from 393873 to 484598. Inverse variance weighted (IVW) method was used as the primary MR analysis methods, and MR Egger, weighted median was performed as supplements. Sensitivity analysis was employed to test for heterogeneity and horizontal pleiotropy. Finally, multivariable MR was utilized to adjust for the influence of BMI on OA. Results After conducting multiple tests (P<0.0023) and adjusting for BMI, MR analysis indicated that a lower AFB will increase the risk of OOA (odds ratio [OR] = 0.97, 95% confidence interval [CI]: 0.95–0.99, P = 3.39×10−4) and KOA (OR = 0.60, 95% CI: 0.47–0.78, P = 1.07×10−4). ALB (OR = 0.61, 95% CI: 0.45–0.84, P = 2.06×10−3) and Age AFSI (OR = 0.66, 95% CI: 0.53–0.82, P = 2.42×10−4) were negatively associated with KOA. In addition, our results showed that earlier AMP adversely affected HOA (OR = 1.12, 95% CI: 1.01–1.23, P = 0.033), and earlier ASOC promote the development of OOA (OR = 0.97, 95% CI: 0.95–1.00, P = 0.032) and KOA (OR = 0.58, 95% CI: 0.40–0.84, P = 4.49×10−3). ALB (OR = 0.98, 95% CI: 0.96–1.00, P = 0.030) and AFSI (OR = 0.98, 95% CI: 0.97–0.99, P = 2.66×10−3) also showed a negative association with OOA but they all did not pass multiple tests. The effects of AAM and NLB on OA were insignificant after BMI correction. Conclusion This research Certificates that Early AFB promotes the development of OOA, meanwhile early AFB, ALB, and AFSI are also risk factors of KOA. Reproductive factors, especially those related to birth, may have the greatest impact on KOA. It provides guidance for promoting women’s appropriate age fertility and strengthening perinatal care.


Introduction
Osteoarthritis (OA) is the most common form of degenerative joint disease, affecting approximately 302 million people worldwide [1].It often leads to long-term chronic pain and impaired joint mobility, significantly impacting the quality of life of those affected [2].With an aging population and changing lifestyle, the increasing prevalence of OA, primarily due to age and obesity, has emerged as a significant public health concern [3].
The incidence of OA varies between genders.Epidemiologic studies indicate that OA is significantly more prevalent in females than in males, with the World Health Organization reporting incidence rates of 18% in women compared to 9.8% in men [4,5].Unique physiological characteristics associated with menstruation and reproduction in women suggest a possible link to OA. Notably, the incidence of OA significantly increases after menopause, indicating that changes in menstrual reproductive factors (MRF) may influence the development of OA [5][6][7][8].One prospective cohort study showed that a later age at menarche (AAM) was related to a lower risk of developing OA, and that past users or users of oral contraceptives or taking hormone replacement therapy increased the risk of joint replacement in OA [9].Riyazi's study showed that later age at menopause (AMP) decrease the risk of familial OA [10].However, there are some different views believe that later AAM Increases the likelihood of knee and hip replacement surgery [11], and menopausal status or AMP are not a predictor of developing OA [12].The Number of live births (NLB) has also been studied, with mixed finding on its impact on OA [9,[13][14][15].
Despite numerous observational studies demonstrating a potential association between MRF and OA, the precise nature of this association requires further clarification.Observational studies are subject to confounding effects and horizontal pleiotropy, which can obscure the distinction between reverse causality and lead to variable results, thus undermining their overall credibility.Mendelian randomization (MR) offers a more robust approach.This analytical method examines potential causal associations between dangerous factors and outcomes using data from genome-wide association studies (GWAS) [16].The genetic variants associated with exposure were used as instrumental variables (IVs) to match the outcome, and MR analysis allowed for a more intuitive response to reflect causal relationship [17,18].By adhering to the principles of random assignment of Mendelian gametes and free combination, MR effectively mitigates issues of confounding factors and reverse causality [17].
In this study, we used a two-sample MR method to assess the possible causal relationship between MRF and OA, and adjust for the influence of BMI, with a view to providing guidance for the prevention and treatment of female reproductive health and OA.

Study design
This two-sample MR analysis consists of two steps.In the first step, we validate the causal association between seven exposure factors related to menstruation (AAM, AMP) and reproduction (Age at First Birth,AFB; Age at Last Birth, ALB; NLB; Age First Had Sexual Intercourse,AFSI; Age Started Oral Contraceptive Pill, ASOC) and three outcomes (Overall osteoarthritis, OOA; knee osteoarthritis, KOA; and hip osteoarthritis, HOA) using a univariate MR (UNMR) method.In the second step, we adjust the effect of Body Mass Index (BMI) on OA through a multivariate MR (MVMR) method.We use Single nucleotide polymorphisms (SNPs) that are significantly correlated with MRF as instrumental variables (IVs).Additionally, in step one, sensitivity analysis was employed to test the stability of IVs,including Cochran Q test, MR-Egger regression, funnel plot, and the leave-one-out method.This design approach has been referenced in several articles [19,20].For SNPs to serve effectively as IVs, they must fulfill the following three assumptions: (1) the relevance assumption, where SNPs are strongly correlated with MRF; (2) the independence assumption, where SNPs are not associated with confounders; and (3) the exclusionary assumption, where SNPs are not straightforwardly involved in the outcome but operate through the way of exposure [21].This study design is shown in Fig 1.

Data sources and SNP selection for exposures and outcomes
All the pooled data for this two-sample MR study were sourced from the publicly accessible IEU GWAS database summary website (https://gwas.mrcieu.ac.uk/).The selected exposures and outcomes originated from different datasets to prevent sample overlap.The exposure data were exclusively derived from the UK Biobank database, while the outcome data came from the EBI database.Both databases primarily encompass European populations.A detailed summary of this information is presented in Table 1.
In order to ensure the selected IVs are representative and minimize bias, we adopted the following criteria for selecting SNPs: We screened for SNPs with a strong significant association (P<5.0×10−8 ) with the exposure and removed those in linkage disequilibrium (r 2 = 0.001, kb = 10000).We used the LDTrait tool (https://ldlink.nih.gov/?tab=ldtrait) to complete the independence and exclusivity assumptions [22].After the initial screening of IVs for exposures we matched the SNPs information for the outcomes to obtain SNPs of outcomes.And proxy SNPs were not used.SNPs with palindromic structures and minor allele frequency (MAF) greater than 0.01 were eliminate.We calculated the F-statistic to measure whether the selected SNPs were weak instrumental variable.an F<10, indicates the presence of a weakly instrumented variable.The F-value is calculated from the following formula: [23], where N is the sample size of the exposed database, k is the number of selected SNPs, and R 2 is the proportion of variance explained by each SNP.R 2 is further calculated as: R 2 where EAF is the effect allele frequency, Beta is the allele effect value, and SE is the standard error.palindromic SNPs and incompatible SNPs.The primary analysis technique employed was the Inverse Variance Weighted (IVW) method.It based on the assumption that SNPs are valid IVs.This method examines the causal association by conducting a meta-analysis of each Wald ratio for the selected SNPs [24,25].MR-Egger and Weighted Median (WM) served as complementary methods of IVW method.MR-Egger presumes that each SNP has an uncorrelated association with horizontal pleiotropic effects, which may lead to inaccurate results, especially when the quantity of genetic instruments is limited [26].The WM method, calculates the causal effect using the median of the estimates and remains consistent even if up to 50% of the SNPs are derived from invalid IVs [27].The IVW method has the highest test efficacy of all the analyses, so we use the IVW approach as the main method of analysis, and we consider such results to be significant as long as the IVW results are significant and the direction of the β values of the other two methods is the same.The MR-Egger method has a low test efficacy, but its intercept value can be used as an assessment of horizontal multivariate validity.To adjust for BMI effect in our models, we conducted Multivariate MR (MVMR) analysis, which also used the IVW method [28].It involves selecting all exposures for SNPs and regressing them against the outcome collectively, with weights determined by the inverse variance of the outcome.This method helps limit the indirect pathway effects of exposure-SNP on other risk factors.The TwoSampleMR (version 0.5.7)package in R software (version 4.3.2) was used as the main analytical tool.A P value of less than 0.0023 (0.05/21, with Bonferroni correction) was considered significant.

Sensitivity analysis
After selecting the IVs and deriving the statistical results from the MR analysis, we also need to perform a sensitivity analysis on the results.In this study, Cochran's Q test, MR-Egger regression analysis, funnel plots, and the leave-one-out method were employed for this purpose.Cochran's Q statistics were utilized to explore heterogeneity among IVs [29].If the Cochran's Q assay is shown to be heterogeneous, we transitioned from a fixed effect to a random effects model, as the latter is more tolerant of heterogeneity.Funnel plots served as a visual tool for observing horizontal pleiotropy, with symmetry in the graphs indicating a lower likelihood of pleiotropy [27].The MR-Egger method was also employed to assess horizontal pleiotropy.A significant deviation from zero in the intercept term of the MR-Egger regression signals horizontal pleiotropy [30].Additionally, the leave-one-out method was used to determine whether the statistical results were influenced by a single SNP by sequentially removing each SNP in turn.In these tests all P-value was employed by 0.05 as statistically significant results.

Results of instrumental variables selection
After filtering SNPs by setting the P-value less than 5×10 −8 and removing the interference of linkage disequilibrium, as well as harmonizing the OA information from the GWAS database, we ultimately identified 523 SNPs as IVs (S1 Table ).These included 194 for AAM, 109 for AMP, 17 for AFB, 6 for ALB, 3 for NLB, 190 for AFSI and 4 for ASOC.The minimum F-statistic value was 32.16, and the maximum was 95.39 (as detailed in S1 Table ).All values were greater than 10, indicating that the IVs could strongly predict OA in this study.

Result of sensitivity test
The Cochran Q test showed that IVs of AAM, AMP, AFB, AFSI and ASOC have heterogeneity.Consequently, we opted for the random effects model of the IVW method (Table 2).The funnel plot shows that all graphs are symmetric, suggesting no significant outliers (S1 Fig) .The MR-Egger results demonstrated that all p-values were bigger than 0.05, demonstrating the effect of instrumental variables away from horizontal polytropy and can adequately represent the exposure factors without interference from confounders (Table 2).The leave-one-out method confirmed that after sequentially rule out each SNP, the estimates remained consistent with the value calculated using all SNPs, proving that the IVs were robust with little influence from individual SNPs (S2 Fig) .

Discussion
In this study, we employed two-sample MR methods to investigate the causal association between two menstrual factors, five reproductive factors and three types of OA.Recognizing that BMI is the most common confounder of OA, we adjusted for its effects using a multivariate MR method.After conducting multiple tests, three exposure factors were conclusively identified as having a significant causal association with OA: AFB, ALB and AFSI.
Both AFB and ALB are negatively causally associated with KOA, and additionally, AFB is also negatively associated with OOA.Although AFB and ALB represent different aspects of reproductive age, a smaller ALB also implies an earlier age of childbearing for the same number of children, thereby giving both exposure factors equivalent meanings.This suggests that an earlier reproductive age increases the risk of both KOA and OOA.This aligns with the findings of some observational studies.For instance, a study involving 28 non-pregnant female participants found that the parous group exhibited smaller knee and hip moments and greater knee flexion angles compared to the nulliparous group [31].This altered habitual loading pattern increases irritation of the cartilage, thereby heightening the risk of developing OA.This suggests that during normal delivery, the expulsion of the fetus affects the vaginal and pelvic structures of the mother, which may, in turn, lead to changes in the force lines of the lower extremities.The alteration in the force lines of the lower limbs may elevate the risk of OA [32,33].Early childbirth could potentially lead to early knee joint injuries, and a greater number of births may lead to significant changes in the force lines of the lower limbs.Another cohort study of 4.6 million Danish people indicated that the risk of OOA in women increased by an average of 1.05 times for each additional child, with 1.10 times increase in the risk of KOA [15].A longitudinal observational study of 1,618 OA patients aged 50-79 years demonstrated that the peril of OA increased with the amounts of births, particularly in women who had more than four children [34].While our results show great similarity with these studies, the NLB was not significant after BMI correction and multiple tests.However, the genetic estimation results of NLB indicated a positive association with OA (OR = 1.03,P = 0.023).Regarding the effect of parity on OA, Wei offered an alternate explanation by suggesting that an increase in number of productions is associated with a decrease in the amount of cartilage in various knee compartments and an increase in cartilage defects [35].Sexual Intercourse has been linked to pregnancy and childbirth, with earlier sexual activity potentially indicating an earlier pregnancy or birth.A study involving 10164 college students show that early sexual behavior promotes the risk of unintended pregnancy, with an astonishing 34.03% of participants experiencing unintended pregnancies [36].While there are no direct observational studies establishing that early sexual behavior increases the risk of OA, related research has shown that early sexual intercourse can exacerbate depressive behaviors [37].This is significant as internalized behaviors like adolescent anxiety, depression, and withdrawnness are known to promote early sexual activity [38].Early sexual debut heightens the risk of teenage pregnancies, births, and miscarriages [39], which indirectly become fertility factors influencing earlier age and an increased number of births.This, in turn, affects the prevalence of OA.
Many experts argue that the decline in postmenopausal sex hormone levels plays a pivotal role in the formation of OA [40,41].They suggest that postmenopausal hormone replacement therapy can mitigate the progression of OA [42,43].However, earlier use of hormone therapy or oral contraceptives does not necessarily confer protection against OA.Observational studies have reported an elevated risk of total knee replacement among individuals who had previously used oral contraceptives or undergone hormone replacement therapy [9].Our findings similarly indicate that earlier ASOC use may tend to accelerate the development of both KOA and OOA.This earlier use of hormone therapy or ASOC might reflect an earlier decline in sex hormone levels in the body.
The impact of AAM on KOA ceased to be significant after adjusting for BMI, indicating a strong correlation between earlier AAM and obesity.This relationship has been substantiated by prior studies demonstrating that higher body weight increases the risk of early menarche in women [44,45].Owing to the inability of observational studies to fully exclude the influence of confounding factors, the conclusion that earlier AAM promotes the risk of OA, especially when considering the combined effects of BMI and AAM, has been drawn from prior research [9,13].Another cross-sectional study supports our hypothesis, confirming that AAM is not associated with the risk of OA after adjusting for confounding factors [14].This suggests that AAM may not be an independent risk factor for OA but may influence OA development as a consequence of obesity.Numerous studies have indicated that the morbidity of OA significantly increases in women post-menopause [46,47].Our study's results show that an earlier AMP does not elevate the risk of OOA, and although it retains a significant association with HOA after BMI correction, it did not remain significant after multiple testing.Nevertheless, our findings provide insight into predicting the relationship between AMP and the development of HOA.
Although only AFB, ALB and AFSI were significantly correlated with OA after adjusting for BMI effects and multiple testing, there is unity in all menstrual and reproductive factors on what they illustrate.These individual factors are closely interconnected, where one may be a prerequisite for another.For instance, individuals with earlier AFSI tend to have earlier ASOC and AFB, and those with an earlier AAM may also engage in sexual activity sooner.There is a strong correlation among AAM, AFSI, AFB and NLB [48].Our results suggest that reproductive factors (AFB, ALB, NLB, AFSI and ASOC) tend to impact OOA and KOA, while menstrual factor (AMP) may influence HOA.The exposure factors we studied (except NLB) almost universally displayed a negative correlation with the outcome.This indicates that later sexual intercourse and later age at childbearing are connecting with a lower danger of OA.These findings are instructive in advocating for avoiding early sexual activity and promoting age-appropriate childbearing.
The advantage of our study lies in utilizing a two-sample MR design and correcting for the influence of BMI using the MVMR method.This approach minimizes the effects of reverse causality and confounding.Additionally, the use of publicly available large sample GWAS datasets for the analysis enhances the accuracy of the results to a certain extent.However, we must acknowledge several potential shortcomings.Firstly, the subjects in our study were primarily of European origin, which may not fully represent the global population.Secondly, menstrual and reproductive factors are intricately linked to changes in human sex hormones, and our study did not account for the effects of various sex hormones on these factors.Furthermore, some exposures yielded a relatively small number of SNPs, potentially affecting the precision of the results.Future research could benefit from exploring larger and more diverse databases.

Conclusions
In conclusion, our study offers some insights into the impact of ARF on the incidence of OA.We found that relatively early AFB, ALB and AFSI are link with a high risk of developing KOA, and later AFB is associated with a decreased risk of OOA.These findings highlight the potential influence of reproductive life events on OA development.However, further studies are necessary to confirm and extend our observations.
Fig 4 presents a comparison of the results from UVMR and MVMR.

Table 2 . Pleiotropy and heterogeneity tests for the causal association between MRF on OA.
Age at menarche; AMP, Age at menopause; AFB, Age at first live birth; ALB, Age at last live birth; NLB, Number of live births; AFSI, Age first had sexual